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FOURTH QUARTER SECOND YEAR PROGRESS REPORT 

This report covers technical progress during the fourth quarter of the second year of 
NASA Sun-Earth Connections Theory Program (SECTP) contract “The Structure and 
Dynamics of the Solar Corona and Inner Heliosphere,” NAS5-99188, between NASA and 
Science Applications International Corporation, and covers the period May 16, 2001 to August 
15, 2001. Under this contract SAIC and the University of California, Irvine (UCI) have 
conducted research into theoretical modeling of active regions, the solar corona, and the inner 
heliosphere, using the MHD model. 

Press Release at the Spring AGU Meeting, Boston, MA, May 30, 2001 

Zoran Mikic was invited by the AAS/SPD press officer. Dr. Criag DeForest, to give a 
presentation at a press conference about the total solar eclipse that was expected in Africa on 
June 21, 2001. An accompanying press release is contained in Appendix A. 

Presentations at the Spring AGU Meeting, Boston, MA, May 29-June 2, 2001 

Our progress during this reporting period is summarized by the following papers that 
were presented by SAIC staff members at the Spring AGU Meeting, Boston, MA, May 29- 
June 2, 2001. 

Predicting the Structure of the Solar Corona During the 21 June 2001 Total Solar Eclipse 
Z. Mikic, J. A. Linker, R. Lionello, and P. Riley (SAIC, San Diego) 

We describe the application of a three-dimensional magnetohydrodynamic (MHD) model 
to the prediction of the structure of the corona during the total solar eclipse that is expected to 
occur on 21 June 2001. The calculation uses the observed photospheric radial magnetic field as 
a boundary condition. This model makes it possible to determine the large-scale structure of 
the magnetic field in the corona, as well as the distribution of the solar wind velocity, plasma 
density, and temperature. We will use magnetic Helds observed on the solar disk prior to 
eclipse day to predict what the corona will look like during the eclipse. The estimated coronal 
density will be used to predict the plane-of-sky polarization brightness prior to the eclipse. 

A copy of this presentation appears in Appendix B. 


Are There Two Classes of Coronal Mass Ejections? A Theoretical Perspective * 

J. A Linker, Z. Mikic, R. Lionello, and P. Riley (SAIC, San Diego) 

Coronal mass ejections (CMEs) are generally accepted as the cause of nonrecurrent 
£ gnetic storms at Earth. Statistical compilations of CME events have shown that CMEs 
launched in the corona can have a wide variation in speeds (Hundhausen et al ., JGR 99, 6543, 
1994). The speed of the CME at Earth and the presence or absence of an interplanetary 
shockwave is an important component of the geoeffectiveness; therefore, the mechanism(s) by 
which fast CMEs might be produced are considered to be of particular importance. Recently, 
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the examination of the acceleration profiles of CMEs has led to the possible classification of 
CMEs as either constant speed CMEs or constant acceleration (St. Cyr et al . , JGR 104, 12493, 

1999; Sheeley et al., JGR 104, 24739, 1999). In this talk, we will examine reasons why one 
might expect CMEs to show two classes of acceleration profiles, and if they do, what the 
implications are for models of CME initiation. Examples from MHD simulations of CMEs for 
different initiation mechanisms will be used to illustrate essential points. 

*Research supported by NASA and Boston University’s Integrated Space Weather Modeling 
project (funded by NSF). 

A copy of this presentation appears in Appendix C. 


Using Global MHD Simulations to Interpret In Situ Observations of CMEs 

P. Riley, J. A. Linker, R. Lionello, Z. Mikic (S AIC, San Diego) 

D. Odstrcil, V. J. Pizzo (SEC, Boulder) 

T. H. Zurbuchen (U. Michigan) 

D. Lario (JHU/APL) 

In this study, we combine two MHD models to simulate the initiation, propagation, and 
dynamic evolution of flux-rope-like CMEs through the corona and out to 1 AU. The coronal 
model encompasses the region of the solar corona from IR S to 20R S , while the heliospheric 
model encompasses 20R S to 1 AU. The CME initiated in the corona propagates smoothly 
across the outer boundary of the coronal solution and through the inner boundary of the 
heliospheric solution. The model solutions show a rich complexity, which, given the relative 
simplicity and idealization of the input conditions, bear a strong resemblance to many observed 
events, and we use the simulation results to infer the global structure of some of these 
observations. In particular, we highlight an event that was observed by both Ulysses and ACE 
in February/March, 1999. At this time, Ulysses was located at - 5 AU and 22°S heliographic 
latitude; thus the two spacecraft were separated significantly both in heliocentric distance and 
latitude. We also use these simulations to separate dynamical effects from force-free models of 
flux ropes in the solar wind. 

A copy of this presentation appears in Appendix D. 


Modeling ofTransequatorial Loops with MH4D 
R. Lionello and D. Schnack (SAIC, San Diego) 

MH4D (Magnetohydrodynamics on a TETRAhedral Domain) is a new algorithm to 
perform magnetohydrodynamic (MHD) simulations of the active regions of the Sun, including 
the large scale coronal structure that surround them. MH4D is a massively-parallel, device- 
independent numerical code for the advancement of the resistive and viscous MHD equations 
on an unstructured grid of tetrahedra. The use of an unstructured grid allows us to increase the 
resolution in the regions of physical interest. A variational formulation of the differential 
operators ensures accuracy and the preservation of the analytical properties of the operators 
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(V-B = 0, self-adjointness of the resistive and viscous operators). The combined semi-implicit 
treatment of the waves and implicit formulation of the diffusive operators can accommodate the 
wide range of time scales present in the solar corona. The capability of mesh refinement and 
coarsening is also included. A preliminary result is presented: a simulation of transequatorial 
loops that include fine details of two interconnected active regions. 

A copy of this presentation appears in Appendix E. 
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APPENDIX A 


Press Release at the Spring AGU Meeting 
Boston, MA, May 30, 2001 
Presented by Zoran Mikic 


4 



May ^y, ^uui 


PREDICTING THE STRUCTURE OF THE SOLAR CORONA 
DURING THE 21 JUNE 2001 TOTAL SOLAR ECLIPSE 

Total solar eclipses offer an excellent opportunity to observe the solar corona. During a total solar 
eclipse the moon blocks the bright light from the solar disk, allowing the faint light scattered by the solar 
corona, which is more than a million times fainter than the photosphere, to become visible. During totality 
the structures that characterize the white-light corona become apparent, including prominences, helmet 
streamers, polar plumes, and coronal holes. Observers who witness a total solar eclipse invariably report 
that it is a beautiful sight to behold. 

On 21 June, 2001 a total eclipse of the Sun will be visible in the southern hemisphere, beginning in 
the South Atlantic, crossing southern Africa and Madagascar, and terminating in the Indian Ocean. Drs. 
Zoran Mikic, Jon Linker, Pete Riley, and Roberto Lionello, of Science Applications International 
Corporation (SAIC) in San Diego, California, have developed a theoretical model to predict what the solar 
corona will look like during forthcoming total solar eclipses. Their model has been used to predict the 
shape of the corona during the eclipse that is expected to occur on 21 June. 

The model is based on the three-dimensional magnetohydrodynamic (MHD) equations that describe 
the interaction of the solar wind with coronal magnetic fields. The group’s results, which have been 
financially supported by the National Aeronautics and Space Administration (NASA) and the National 
Science Foundation (NSF), are being presented by Dr. Mikic at the American Geophysical Union meeting 
in Boston. 

The calculation relies on Earth-based measurements of the magnetic field in the solar photosphere to 
infer the structure of the solar corona. The measurements are taken at the National Solar Observatory at 
Kitt Peak . “It is remarkable that measurements of the magnetic field in the photosphere can tell us so 
much about the corona,” says Dr. Mikic. “The simulations already have a strong resemblance to coronal 
images. In the future, the agreement will only improve as we refine the physics in our model and as we 
utilize ever-faster computers. We will be able to study even finer details in the corona.” The researchers 
use supercomputers at the San Diego Supercomputer Center (part of NSF’s National Partnership for 
Advanced Computing Infrastructure) and the Department of Energy’s National Energy Research 
Supercomputer Center to solve the equations. The output from the model is used to predict the brightness 
of the corona. 

So far, this model has been applied to five eclipses (see http ; / /haven . saic . com) . This kind of 
modeling helps us to understand the structure of the solar corona, especially the location of helmet 
streamers and coronal holes, and the nature of the fast and slow solar wind, and provides a rudimentary 
test of predictive capability. Eclipses offer an opportunity to test such models, and to understand the 
influence of the Sun on the Earth. 



Predicted polarization brightness (pB) for the 21 June, 2001 eclipse, together with traces of the magnetic field lines, at 13: 10 
UT (corresponding to totality in Lusaka, Zambia). The pB signal is produced by white light scattered off electrons in the 
coronal plasma. This is the view of the Sun that would be seen by an observer on Earth with a camera aligned so that 
vertical is toward the Earth’s north pole. 

Zoran Miki6, SAIC, 10260 Campus Point Drive, San Diego, CA 92121; Tel: 858-826-6934; Email: mikicz@saic.com 
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APPENDIX B 


“Predicting the Structure of the Solar Corona During the 21 June 2001 Total Solar Eclipse” 
Z. Mikic, J. A. Linker, R. Lionello, and P. Riley 
Presented at the Spring AGU Meeting 
Boston, MA, May 29-June 2, 2001 
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Introduction 


• The solar magnetic field plays a key 
role in determining coronal 

• The principal input to MHD models is 
the observed solar magnetic field 

• 3D MHD models can be used to 
compare with eclipse and coronagraph 
images, SOHO images (LASCO, EXT), 
Ulysses and WIND spacecraft data, 
and interplanetary scintillation (IPS) 
measurements 

• MHD computations can tell us about 
the structure of the corona 

• Eclipses can help us to verify the 
accuracy of the models 

• 21 June, 2001 total solar eclipse: 
visible in the southern hemisphere, 
(South Atlantic, southern Africa, 
Madagascar, and Indian Ocean) 

• Totality in Lusaka, Zambia is at 13:10UT 
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The Polytropic Model 

• Neglect thermal conduction, coronal 
heating, radiation loss, and Alfven 
waves (set p w = 0 and S ~ 0) 

• Simulate these effects (crudely) by 
setting y = 1.05 (Parker 1963) 

• A possible extension is to have y = y(r), 
with 7 increasing far from the Sun 

• We have used this model extensively in 
3D computations of the structure and 
dynamics (e.g., CMEs) of the solar 
corona 

• The corona is modeled reasonably well, 
but the properties of the interplanetary 
solar wind are not accurate (speed, 
density, temperature) 

• An improved model (with more accurate 
energy transport) is being developed 


MHD Equations 
(Polytropic Model) 
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Finding MHD Solutions 


• Use line-of-sight magnetograms to 
deduce B r from B\ os at r = R Q 

(e.g v Kitt Peak Solar Observatory and 
Wilcox Solar Observatory synoptic 
maps) 

• Calculate a potential field matching B r at 
r = R 0 

• Specify T and p on the solar surface 
r = R 0 (e.g., uniform T 0 and p 0 ) 

• Set up p, p, and v from a spherically 
symmetric solar wind (ID Parker 
solution) 

• Integrate 3D MHD equations in time 
until steady state is reached 

• This gives the structure of the coronal 
magnetic field B (as well as p, p, v, T) 

• Compare with observations 


Polarization Brightness 

• Light scattered off the coronal electrons 
is observed in coronagraphs 

pB(x) - K J n e (x-x')C(r') dl' 

los 

• C(r) is a scattering function (e.g., Billings 
1966) 

• To produce a plane-of-sky image, we 
apply a (radial) filter to pn ("vignetting 
function") and we simulate the effect of 
an occulting disk 


"Canonical" Helmet Streamer 

• Idealized helmet streamer configuration 

• Start with 2D (axisymmetric) dipole field 
(Pneuman & Kopp 1971) 

• Closed-field region with a static (v = 0) 
dense plasma, surrounded by an open- 
field region with solar wind streaming 
along the field lines 

• A current sheet surrounds the helmet, 
and an equatorial current-sheet 
separates fields of opposite polarity 


Whole Sun Month 

10 - Sep. 8, 1996 


Radial Velocity 
Open and Closed Field Lines 
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Comparison with Eclipse and 

CORONAGRAPH OBSERVATIONS 


• November 3, 1994 eclipse (Chile) 

Compare with eclipse image (HAO) and 
Mauna Loa coronagraph data 

• October 24, 1995 eclipse (Vietnam): 

We predicted the coronal structure on 
Oct. 5, 1995 

Compare with eclipse image 
(S. Koutchmy) 

• March 9, 1997 eclipse (Russia, China, & 
Mongolia): 

We predicted the coronal structure on 
Mar. 3, 1997 

Compare with eclipse image (E. Hiei) 

• February 26, 1998 eclipse (Carribean): 

We predicted the coronal structure on 
Feb. 13, 1998 

Compare with eclipse image (HAO) 





August 11, 1999 eclipse (Central 
Europe, Turkey, Iran): 

We predicted the coronal structure on 
July 28, 1999 

Compare with eclipse image (Espenak) 










Comparison of a 3D MHD Coronal Prediction with 
an Image of the 11 August 1999 Total Solar Eclipse 


Evolution of the Photospheric Magnetic Field Leading 
up to the August 11, 1999 Total Solar Eclipse 

(Courtesy of Kitt Peak Solar Observatory) 

CR1949 (May 1 - May 28, 1999) 

tjT T p »* J 


Fred Espenak’s Composite Image (Turkey) 


Longitude 


CR1950 (May 28 - June 24, 1999) 


Ap; 


Predicted Polarization Brightness (MHD Model) 




Predicted Magnetic Field Lines (MHD Model) 

Figure 1. Comparison between a composite eclipse image created from photographs taken by Fred 
Espenak in Lake Hazar, Turkey (top) with the predicted polarization brightness of the simulated 
solar corona from our 3D MHD model (middle). The projected magnetic field lines from the 
model are also shown (bottom). Terrestrial (geocentric) north is vertically upward. The eclipse 
image is copyrighted 1999 by Fred Espenak. 












Latitude (deg.) Latitude (deg.) 


August 11, 1999 Total Solar Eclipse 

Kitt Peak Synoptic Chart (CR1951), B r 



Smoothed Magnetic Field (used in MHD model) 



Carrington Longitude (deg.) 


June 21, 2001 Solar Eclipse 
Prediction 

• Started a high-resolution simulation on 
May 19, 2001 using magnetic field data 
from a Kitt Peak Solar Observatory 
synoptic map of CR1975 (April 9-May 
6 , 2001 ). 

• Low-resolution case: 61 x 71 x 64 
(r,0,0) mesh points 

(~ 7 CPU hrs on the Cray-T90) 

• High-resolution case: 111 x 101 x 128 
(r,0,0) mesh points 

(~ 90 CPU hrs on the Cray-T90) 

• This calculation will be updated in the 
next few weeks with more recent 
magnetic field data 

• The code runs on PCs and Cray 
computers. It is also being developed 
to run on massively parallel machines 
using MPI (e.g., Cray T3E, Beowulf, IBM 
SP3) 

• See: 

http: / / haven.saic.com / corona / modeIing.html 



Predicting the Structure of the Solar Corona 
During the 21 June, 2001 Total Solar Eclipse 
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Predicted polanzntion bnghtness (lop left) together with traces of the magnetic field linn id the solar coma (top right) fee the eclipse expected on 21 June. 3001 
st 13 10 UT (with terminal north up) The Sun's surface shorn s color contours of the radial component of the measured photospbenc magnetic Geld [root Kill 
tValr National Sri nr nt wervatorv. storing lltc location cf active reports (strong magnetic fields) Click the images for higher resolution pictures To new these 
images in a coordinate system aligned with solar north, dick here: 
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'Hie photosphenc magnetic field maps « e use for our calculations are bull up from 
dails observations of die Sun dun rig a solar rotation I he sc maps give a good 
approximation of the Sun's magnetic flux if the large-scale flux is not changing much 
throughout a rotation Res tously. w e have computed coronal models for an eclipse 
during the declining phase of the last solar cy cle t November 3. 19M> . far In o eclipses 
during solar minimum f Ocloher 24. 1993 and March 9. 1997) - one eclipse diking the 
the early rising phase of solar cycle 23 ( February 26. 19981 . and one edipsc 
approaching solar maximum ( August I I. 19991 The June 21 ,2001 eclipse, ntuch is 
occurring at the peak of solar maximum, presents a challenge. The photosphenc 
magnetic field is eioiving more rapid!}, making synoptic magnetic field data a less 
reliable approximation to the line slate of the photos poetic magnetic field The 
complexity of coronal structures requires high resolution rum (requiring approximately 
1 .500,000 gnd points in our calculation) 

These figures show the photosphenc magnetic field maps for three Carrington rotations, 
CR 1975, CR 1976, andCR197f>+1977, as measured by the National Solar Observatory 
at Kin Peak . The maps show the measured photosphenc magnetic field ax a function of 
latitude (vertical axis) and Carrington longitude (horizontal axis) Red show s outward 
directed magnetic flux, and blue shows inward directed fiux dick the images for 
higher resolution pictures These maps are considerably more complex than mapa 

daring niat minimum 


CRI975 (April 9 • May 6, 2001) 



s. i 

CR1976 (May « - June 2, 2001) 



CR1974+1977 (May 19 - Juna 15, 2001) 



Movies: 

We hare made a movie of the pdan ration brightness from our MHD simulation of the solar corona during Canington rotation 1975 (April 9 ■ May 6, 2001) This 
illustrates visually how rapidly the solar corona changes as a result of solar refahoo during the maximum phase cf the solar cycle Ycrn can get an MPFn vrrxim 
(599 kbytes) or a QuickTime version ( 1 .003 kbytes). 

If your movie player can continuously loops movie while playing it ael Ihn option to ■on’ for the tea effect Foresample on a Silicon Graphics machine, you would use the execute line 
m 03 •sepia) er A 0 cr 1975 mpf 


Other web resources for the eclipse: 

• Fred Esnenak's 2001 Total Solar Felipse website (A meat resource for the eclipse, including lots of information and links to other sites) 

• Eclipse information from N ASA’s Goddard Space Flight Center 

• Star and TdncaxMaawnc 

• 1AU Waking Group on felipsn 

• MrEdmse.com 

• Upcoming Edumts 


Return to the Coronal Modeling Page 
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APPENDIX C 

“Are There Two Classes of Coronal Mass Ejections? A Theoretical Perspective” 
J. A Linker, Z. Mikic, R. Lionello, and P. Riley 
Presented at the Spring AGU Meeting 
Boston, MA, May 29-June 2, 2001 



ARE THERE TWO 
CLASSr s OF CORONAL 
MAS 3 EJECTIONS? 

A THEORETICAL 
PERSPECTIVE* 


JON A. LINKER 
ZORAN MIKIC 
ROBERTO LIONELLO 
PETE RILEY 

Science Applications International Corporation 


*Research Supported by NASA and NSF (through 
Boston University's Integrated Space Weather Modeling 
Project). Computations performed at NPACI/SDSC. 


Introduction 


• Despite many years of study, the origin 
evolution and of coronal mass ejections 
(CMEs) is poorly understood. 

• What are the underlying physical 
differences between "fast" and "slow" 
CMES? => Important for Space Weather 

• Recently, it has been demonstrated that 
many CMEs observed with the LASCO 
coronagraph fall roughly into two classes 
(Sheeley et al. 1999): 

(1) CMEs that accelerate up to an 
asymptotic speed; 

(2) CMEs that travel at constant speed 
or decelerate 



TimeHeight Plot (LASCO DATA) 


Ni-PQS 



Sheeley et al., 1999 


From Sheeley et al., 1999 
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My "Assignment:" 


"The modelers should address the question 
of whether there is any theoretical reason 
why an accelerating CME requires a 
different expulsion mechanism or 
environment than a CME that starts out fast 
and then decelerates. 

They should ask themselves whether they 
can produce the different speed profiles just 
by varying the source region. ... 

Of course they should also suggest other 
reasons they can think of that might result 
in the different speed profiles. " 
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Questions 

• Does the division of CMEs into 
"accelerating" and "constant speed" 
imply more man one mechanism for 
CMEs initiation? Does it imply only one 
mechanism? 

• Can we use this classification to help us 
to understand CME initiation or 
constrain the possible mechanisms? 


Mechanisms for Solar 
Activity 

• CMEs, fkr;^, and prominence eruptions 
require si^Jficant amounts of energy 
(~10 32 ergs) 

• Most theories assume that energy is 
released from the coronal magnetic field 
(see next talk for an opposing view) 

• There are many observations of 
nonpotential magnetic field structures in 
the corona harboring significant amounts 
of magnetic energy 

• What causes this energy to be released? 
Many candidate mechanisms 

• What sort of time-height profiles of 
ejecta are implied by different 
mechanisms r Today we consider two: 

(1) Shearing of the photospheric 
magnetic field 

(2) Magnetic flux cancellation 


Eruption by Photospheric 
Shearing Flows 


• Previously, we showed that photospheric 
shear leads to the eruption of magnetic 
field arcades and helmet streamers if the 

f >hotospheric shear exceeds a threshold 
Mikic &Linker 1994,Linker &Mikicl995) 

• Start from a helmet streamer 
configuration 

• Introduce flows (typically 0.5-5 km/ s) at 
the photosphere that twist or shear the 
magnetic field and energize the 
configuration. 

• When the magnetic shear crosses a 
threshold, eruption occurs. Eruption 
threshold does not depend on how fast 
the shear is introduced. 

• Eruption is related to Magnetic 
Nonequilibrium - the appearance of a 
discontinuity in force-free equilibrium 
configurations 



i : id Line Evolution 



A^max ~ 1 .60 R$ Ay max = 1 .71 R$ 


_ il 

Eruption of a 3D Helmet Streamer triggered Shearing Flows j j 


Eruption by Flux Cancellation 


• We have found that reduction of 
magnetic flux near the neutral line of a 
sheared or twisted magnetic 
configuration (i.e., flux cancellation) can 
lead to the formation of stable flux rope 
configurations (Amari et al , 1999, 2000; 
Linker et al, 2001) 

• The flux ropes are capable of supporting 
prominence material in the corona 

• When the flux cancellation crosses a 
threshold, the entire configuration 
erupts. In the case of a helmet streamer 
configuration, a CME is ejected into the 
simulated solar wind. 

• The eruption is more energetic than 
eruptions triggered by photospheric 
shearing flows. 


Consider 5 Model Problems: 


Vary mechanism: Shearing and Flux 
cancellation 

Vary initial corona: Base temperature 
1.4X10 6 K or 1.8X10 6 K (solar wind 
speeds of 250 km/s or 350 km/s at 20Rs 

Performed one flux cancellation 
simulation with broader shear 
distribution 


MHD Equations 



Helmet Streamer Configuration for CME Studies 


Flux H'foz) 



Initial Potential Field 


Relaxed Helmet Streamer 


Current Density 




Relaxed Helmet Streamer 





Eruption of a Helmet Streamer 
By Emerging Flux 


Flux'Ffr.z) 



Unsheared streamer Sheared streamer 4i% emerged flux 

t as / 0 t = t Q + 6 hours 



7.5% emerged flux 10.5% emerged flux 12% emerged flux 

t = f D + 10 hours t = t Q + 14 hours t = t 0 + 16 hours 



13.5% emerged flux 
t - 1 0 + 18 hours 
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15 % emerged flux 
t = t Q + 20 hours 


15% emerged flux 
t = t Q + 2.5 days 


Eruption of a Helmet Streamer 
By Emerging Flux 


Polarization Brightness 



Unsheared streamer 

Sheared streamer 

4.5% emerged flux 



t = t Q + 6 hours 



7.5% emerged flux 10.5% emerged flux 12% emerged flux 

t - # 0 + 10 hours t = t 0 + 14 hours t - 1 0 + 16 hours 



13.5% emerged flux 15% emerged flux 15% emerged flux 

t = t 0 + 18 hours / = f 0 + 20 hours / = / Q + 2.5 days 









Latitude Height (R s ) j j Latitude Height (R g ) 


TimeHeight Profiles: Flux Cancellation and Shearing Simulations 


Flux Cancellation 1 



Time (Hours) 


Shearing Motions 1 



TimeHeight Profiles: Slower Solar Wind 


Flux Cancellation 2 Shearing Motions 2 
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TimeHeight Profiles: 


Flux Cancellation 3 Shearing Motions 1 
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Flux Cancellation, Slow Wind Background 


Height vs Time 
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Speed vs. Height 



Shearing Flows, Medium Wind 

Height vs. Time 
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Speed vs. Height 
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Shearing Flows, Slow Wind 
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Summary 


• We examined the results from simulated 
CMEs initiated by magnetic flux 
cancellation and by pnotospheric 
shearing flows. 

• The flux cancellation CMEs yielded 
either accelerating or constant speed 
eiecta depending on the properties of 
tne ambient solar wind. 

• The CMEs initiated by photospheric 
shearing flows yielded accelerating 
ejecta. Shearing flow simulations with 
other initial configurations might also 
yield constant speed ejecta. 


What Have We Learned? 


Questions I raised at the beginning of the 
talk: 

• Does the rough division of CMEs into 

"accelerating" and "constant speed" 
imply more than one mechanism for 
CMEs initiation? No. 

• Does this division imply that CMEs are 
initiated by only one mechanism? 

No. 

• Can this classification by itself help us to 
understand CME initiation or constrain 
the possible mechanisms? No. 

• Our results suggest that the approximate 
division of CMEs into "accelerating" or 
"constant speed" is a natural 
consequence of having a range of energy 
inputs and diverse solar wind conditions 
for different CME ejecta. 
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Overview 


■ Introduction 

• Coronal observations of CMEs 

• In situ observations of magnetic clouds 

• Modeling Flux ropes 

> MHD modelling of CMEs 

• Coronal model (SAIC) 

• Heliospheric model (NOAA/SEC) 

■ Comparison with observations 

• ACE/Ulysses spring 1999 event 

■ Summary 

■ Future work 





Introduction: White Light Observations 

of CMEs 



Introduction: In situ observations of 

magnetic clouds 



Introduction: Models of Magnetic 
clouds/flux ropes 


* Kinematic Models: 

- Force-free 

> Force-free with expansion 
< None force-free with expansion 
Dynamic Models: 

■ Fluid 

- MHD 


MHD Simulation: Coronal solution 
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MHD Simulation: Coronal solution 



MHD Simulation: Coronal solution 
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MHD Simulation: Heliospheric solution 
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MHD Simulation: Heliospheric solution 
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Comparison of an observed magnetic cloud 
with the generic simulation 
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ACE/Ulysses Feb 1999 CME Event 



ACE/Ulysses Feb 1999 CME Event 





Ulysses/ACE May 1999 CME 



Ulysses/ACE May 1999 CME 





Force-Free Fitting at ACE and Ulysses 


Parameter 


Start (day honin' 


Stop (day hr:min 


<l> (deg) 


e 


;/R 


He 


R 


ACE 

Ulysses 

49 14:00 

62 20:00 

50 11:00 

64 22:00 

282.1 

271.1 

-1.3 

53 

0.738 

0.0064 

-1 

-1 

0.28 

0.34 


ACE/Ulysses Feb 1999 CME Event: 
Global structure from force-free fitting 

























Comparison of in situ observations 
with MHD simulation: Ulysses 



Comparison of in situ observations 
with MHD simulation: ACE 





Deceleration of CME 


■ Estimate time for CME to travel from ACE to Ulysses 
Fcme(ACE) = 590 km/s 

^CME(lJlysses) = 46o km/s 
Assume constant deceleration between 1 and 5 AU 
=> At = 13.62 days 

■ CME actually arrived 13.9 days later 

=> Large-scale structure not significantly distorted 


ACE/Ulysses Feb 1999 CME Event: 
Global structure 






Summary 


■ Idealized MHD simulations: 

• provide global context for interpreting in situ CME 
observations 

• Emphasize significant role of dynamics on the 
evolution of ejecta and their associated disturbances 

• Identify new effects that can be searched for in 
observations 


I 

| 

i 
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Future Studies 


■ Explore effects of different launch profiles 

• Do different mec hanisms of eruptions lead to different 
interplanetary signatures? 

■ Launch CME in three dimensions 

• 3-D flow . 

• Magnetic connectivity to Sun 

■ Incorporate composition into model 
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GOALS AND CHALLENGES 


• Develop an efficient 3-D discrete representation 
for the resistive MHD model using an 
unstructured grid of tetrahedral cells 

• Truly arbitrary geometry 

• Use Cartesian coordinates 

• Avoids coordinate singularities and complicated metrics 

• Apply to a variety of problems 

• Fusion 

• Solar and space physics 

• ?? 

• Challenges 

• Discrete representation of differential operators 

• Compactness (couple only nearest neighbor grid points) 

• Self-adjointness 

• Annihilation properties (eg., V B = 0) 

• Solution of implicit system 

• Grid generation 

• Implementation 

• Code and data structure (OO techniques) 

• Parallelism 

• Grid decomposition 


THE FIRST CHALLENGE 

Find efficiently computable discrete 
representations of the magnetic field and 
differential operators on a 3-D unstructured grid 
of tetrahedra 

• The magnetic field and current density: 

B = V x A, J = V x B t 

J = Vx Vx A 

and consequently 

• Both B and J are solenoidai 

• The operator is self-adjoint: 

provided A, and A 2 satisfy the same boundary 

J A 2 • V x V x A ^dV = J A 1 • V x V x A 2 dV 
conditions 

• Seek a discrete representation that retains 
these properties 


GEOMETRY 



- Sides are labeled by the index of their opposite vertex. 

- C is the centroid of the tetrahedron. 

- rrijj is the midpoint of edge i;j. 

- C| is the centroid of side i. 

- Sj is the vector area of side i. 

- mjj C|, C, and are coplanar. 

- s; is the area of the dual median surface. 


Sj = Sj/3 




FINITE VOLUME METHOD 


/ 


• Use integral relations to define differential 
operators 

• Examples: 

Gradient: 


jVfdV fdS 

Divergence: 

Jv-Fc/y = Jji-FdS 

Curl (2-D): 

J^-VxFdS-^F-d 

Curl (3-D): 

xFdV = jpxfdS 

• Integrate over control volume (computational 
cell), eg., 

V C S 

• We will apply this technique to tetrahedral and 
dual median volume elements 


THE MAGNETIC FIELD 


• The vector potential A is defined at vertices of 
tetrahedra 

A varies linearly within a tetrahedron 

• Integral definition: 

J BdV = J V x MV = J^JS x A 

• Apply to tetrahedral cell: 

B t V r = £s s xA s 

S=1 

B r Is constant within tetrahedron z 
A s is the average of A over the 3 vertices of side S 
V T is the volume of tetrahderon z 

• In terms of vertices 

B r=-^-X S v xA v 
dV r v ( t ) 

v(t) are the indices of the 4 vertices of tetrahedron z 
S v is the (outward directed) side opposite vertex v 


DIVERGENCE OF B 


• Apply Gauss’ theorem to dual median volume 
element surrounding vertex v 

v v = X t(v ±1/ t is the volume of this cell 
t(v) are the indices of the tetrahedra sharing vertex v 

= |l4l rW Sr(S,xS r ) 

v(t) is the index of the side of tetrahedron x opposite vertex v 
y( t) are the indices of the 4 sides of tetrahedron x 

• This gives: 

• NOW A MIRACLE OCCURSI!! 

After some algebra , we find that the contributions 
from common faces of adjoining tetrahedra cancel 
exactly! 

• The result! 

(VB) v = 0 

• Caveat This holds for interior vertices only 


ALTERNATE DERIVATION OF B 





THE CURRENT DENSITY 

• Use generalized Stokes’ formula as applied to 
dual median volume element 

\*dV = J V xBdV = jji xBdS 

d r(v) 

r(v)v'(r) 

M T (v, V ) = gyr[(Sv(T) • S^( T ))l - S V -( T )S V ( T )] 

• Substituting the expression for B in terms of A: 

• Symmetry: 

M ap{v,v') = M a,p~x,y,z 

• This is the discrete “curl-curl” operator 


SELF-ADJOINTNESS 


• Define inner product: 

(p,a)= Jp-odv^ -Q>. 

T=1V{T) 4 

• L is self-adjoint if : 

(P,LQ) = {Q,LP) 

• Direct calculation using symmetry of M : 

(P,M Q) = 

= l4^XX-°v'">>)P, 

= (Q ,M P) 

• The discrete operator is self-adjoint 



VARIATIONAL PRINCIPLE 


• Minimize functional 

/(A) = ||[(VxA) 2 -2J-A]cN/ 

• Variation of A: let 

A -> A + e8A 

/(A + eSA) = 1 J [(' V x A) 2 - 2 J . A]d V 

+ eJ(VxAVx<5A-J-5A)cM/ 

+ |e 2 J(Vx5A) 2 dV 

• For minimum, coefficient of e must vanish: 

JSA • (V x V x A - J)dV = -J s <5Af • (B x n)dS 
A t =(l-nn)- A 

• Natural boundary condition: 

SA t = 0, or A t specified on boundary (Dirichlet) 

• Solutions v x v x A = J, with J specified in V, and 
A, specified on the boundaries, minimize < A). 


DISCRETE MINIMIZATION 

• Minimize 1(A) on a grid with N v vertices and N T 
tetrahedra 

• Expand A(x) in basis functions 

A (x) = X A v«v(x) 

v 

• Expand J in delta-functions 

J(x) = lMv5(X-X v ) 

v 

• Substitute into 1(A): 

/ = ^IX[ A V M (v,' v ') • A v , - V v j v ■ A v <ov(x v )] 

V v' 

M ( v, v') = J [(' V a v • V a V ’ )l - V a V ’ S7 a v ]dV 

• To minimize, set 


= 0, y = x,y,z 


• Result: 


T (p,v)-A v , 


• With tent expansion functions gives 

finite volume expression 



BOUNDARY CONDITIONS 

No reference to boundary conditions in discrete 
minimization 

Discrete expression for the “curl-curl” operator 
is 3N V equations in 3N V unknowns 
Could be solved for all unknowns , Including all 
values at the M v boundary vertices 

Absence of surface term implies that solution 
will satisfy the natural boundary condition 

SA r (Bxn) = 0 

Since A f is nof fixed, this can be satisfied only if 

Bxn = 0 

Constraint on source and surface field: 

Solutions with vanishing tangential magnetic field 

exist only If total current vanishes 
In general, we must specify A r on the boundary 


ALTERNATE BOUNDARY 
CONDITIONS 

• Must include surface term in the functional to 
be minimized 

/(A) = J[(V x A) 2 - 2J • A ]cN/ + J(A r P - 2Q) • A t dS 

S 

Pisa self - adjoint matrix and Q is a vector 

• Performing variation and minimization as before 
leads to the condition 

JfiA-(Vx Vx A- J)c/y + JSAf (Bxn+P- A f -Q)c/S = 0 

• Since the tangential variation no longer 
vanishes at the surface, this can be satisfied 
only if 

• Mixed (von Neumann/Dirichlet) boundary 

VxVxA = J/n I/, andBxh+P ■ A t =Q on S 

condition 

• Allows specification of tangential magnetic field 

• Could apply discrete minimization to this 
modified functional to obtain more formalism 
that accommodates these boundary conditions 


MHD FORMULATION 



• Vertices: 

A, J, andpv 

• Centroids: 

p, p, and B 

• Velocity averaged to faces or centroids, as 
required 

• Ail quantities advanced in time by applying 
conservation laws to control volume 

• Use the NIMROD anisotropic semi-implicit 
operator 


SUMMARY AND STATUS 

• Developed a formalism for defining A, B, and J 
on an unstructured grid of tetrahedra 

• Finite-volume approach 

• B is solenoidal 

• “curl-curl” operator is compact and self-adjoint 

• Solutions of discrete equations minimize the 
same functional as solutions of differential 
equations 

• Tangential A can be specified on boundary 

• Can be generalized to specify tangential B 

• Caveats: 

• Demonstrated only first order accuracy 

• B Is solenoidal only on interior vertices 

• Have not demonstrated that J is solenoidal 

• Solved model implicit resistive diffusion 
problem 

• Starting point for full MHD model 



PRELIMINARY VALIDATION 

• Successive application of V x operators on 
primary and dual grids is identical to the 
composite V x V x operator 

i 

• Verified that V • B = 0 

• Solution of diffusion equation in cubic, 
cylindrical, and spherical domains 

• Comparison with analytic solutions 

• Single and multiple processor calculations yield 
identical results 

i 

• Future directions 

• More complex geometry 

• Magnetostatic problem 

• Must deal with gauge condition 

• Reconstruction of force-free fields from boundary data 

• Resistive MHD 


Implementation 


• Standard Fortran 90 

- Object-oriented features 
http://www.cs.rpi.edu/~szymansk/oof90.html 

• Parallelization 

— Runs on serial and parallel machines 

- MPI 

http://www.mpi-forum.org 

• Grid generation with LaGriT (LANL) 
http://www.tl2.lanl.gov/~lagrit 



Mesh Partitioning between 5 PEs 


Implementation (cont.) 


• Grid decomposition with Metis (U.Minn.) 
http://www.cs.unm.edu/~metis 


• Matrix inversion with PETSc (ANL) 
http://www.mcs.anl.gov/petsc 


• Graphics with GMV (LANL) 
http://www-xdiv.lanl.gov/XCM/gmv/ 


i 






Test Case: 
Diffusion Equation 




Numerical vs. Analytical 
Solution 

of Diffusion Equation: Box 

• A solution of the diffusion equation in a 
box domain: 

A z = exp(-i/t) sin sin ^ 


Then 


Lx Ly Tlx fly — 1 , 

1 ) = 0.01 

^Analytical = 0.197. 


• Using 125 nodes, 383 cells in a 1 x 1 x 1 
cubic domain and A t = 0.01, 


^Numerical — 0.182. 





''Numerical — 0.188 


Numerical vs. Analytical 
Solution 

of Diffusion Equation: Cylinder 

• A solution of the diffusion equation in a 
cylinder: 

A z = exp(-i/t)Jo (fcr) 

• If 

k = 2.405 
7] = 0.04 

Then 

^Analytical ~ 0-231 

• Using 500 nodes, 1704 cells in a cylindrical 
domain of h = R = 1 and At = 0.01, 


^Numerical — 0.218 
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Numerical vs. Analytical 
Solution 


of Diffusion Equation: Plots 
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Numerical vs. Analytical 
Solution 

of Diffusion Equation: 
Convergence 



Grid Size 


1000 



Density Advection 


• Advection equation: 


f = -v 

dt 


pv. 


Finite volume formulation with upwinding: 


Ap^ 

At 

nUd) 


Vi 


Es l-iifF®. 


= < 


p 1 if S( ■ Vj > 0 
p J if Sf • Vj < 0 





Momentum Advection 



N0I133AGV WfllN3l/\IOW -QfrHIfll 




Potential Models 


• Active region during Whole Sun Month (August- j 

September 1996). j 

• Active regions AR8102 and AR8100 (Novem- 
ber 1997). 

• We use synoptic magnetograms from Kitt 
Peak National Observatory. 

i 

i 

I 

• The grid extends from 1 to 10 #©. 


• The meshes consist of 42001 (32881) ver- 
tices and 242705 (189237) tetrahedra. 

• The resolution varies from 8 degrees (out- 
side the active regions) to 1 degree (in- 
side). 






High Resolution and Large Scale Model of the Corona (IV) 
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